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Abstract 

A large-eddy simulation (LES) code that utilizes high-resolution numerical schemes is 
described and applied to a compressible jet flow. The code is written in a general manner 
such that the accuracy/ resolution of the simulation can be selected by the user. Time dis- 
cretization is performed using a family of low-dispersion Runge-Kutta schemes, selectable 
from first- to fourth-order. Spatial discretization is performed using central differencing 
schemes. Both standard schemes, second- to twelfth-order (3 to 13 point stencils) and 
Dispersion Relation Preserving schemes from 7 to 13 point stencils are available. The code 
is written in Fortran 90 and uses hybrid MPI/OpenMP parallelization. The code is applied 
to the simulation of a Mach 0.9 jet flow. Four-stage third-order Runge-Kutta time stepping 
and the 13 point DRP spatial discretization scheme of Bogey and Bailly are used. The high 
resolution numerics used allows for the use of relatively sparse grids. Three levels of grid 
resolution are examined, 3.5, 6.5 and 9.2 million points. Mean flow, first-order turbulent 
statistics and turbulent spectra are reported. Good agreement with experimental data for 
mean flow and first-order turbulent statistics is shown. 
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Nomenclature 

finite difference stencil coefficient 

filter stencil coefficient 

shear layer thickness 

phase speed 

numerical phase speed 

total energy 

arbitrary function 

grid indices 

turbulent kinetic energy 

pressure 

heat flux vector 

radial distance from centerline 

time 

velocity tensor 

axial, radial and azimuthal velocity components 
axial, radial and azimuthal turbulence intensities 
coordinate tensor 
cartesian coordinates 
jet diameter 

spatial discretization operator 
filter damping function 
energy spectral density 

Mach number, number of Runge-Kutta stages 
stencil half width 
jet Reynolds number 
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T 

temperature 

Uj 

ideal jet velocity 

At 

time step 

Ax 

grid spacing 

a 

angle of attack 

& rn •> fim 

Runge Kutta scheme coefficients 

K 

wave number 

P 

viscosity 

P 

density 

a 

filter strength coefficient 

Tij 

stress tensor 

subscripts 

jet 

jet exit condition 

0 

stagnation/total condition 

oo 

freestream condition 

superscripts 

- 

spatial filtered 


favre filtered 


filtered 

sgs 

sub-grid scale 


I. Introduction 

The prediction of the flowfield of turbulent jets using computational fluid dynamics (CFD) is an important 
research area. Turbulent jet flows are a critical component in propulsion, aeroacoustics, combustion and 
low-observable (reducing the infrared signature of the exhaust) technologies. A better understanding of the 
flow physics and improved predictive tools could lead to quieter, more efficient and less vulnerable aircraft. 
The standard practice in the aerospace community for predicting these flows is using Reynolds Averaged 
Navier-Stokes (RANS) based CFD codes. RANS methods are based on the time-averaged Navier-Stokes 
equations and rely on a model to represent the effects of the turbulence. Since jet flows are dominated by 
turbulent mixing, the turbulence model in the RANS solver is the critical component in the simulation. 1 
Despite decades of development, to date RANS has failed to produce reliable predictions of the jet flowfield. 2 
Standard models such as Menter’s SST model 3 and k — e 4 fail to accurately predict the potential core length 
(the inviscid portion of the jet not effected by the shear layer) and shear layer growth rate. Turbulence 
models specially formulated for round jets produce better results, but these models can not be generalized 
for other types of jets. 5 In addition, because RANS produces a steady /time- averaged solution, no detailed 
information on the turbulent structures is generated. 

Large-eddy simulation (LES) is a CFD technique that offers the promise of improved predictions and 
increased information for turbulent flows. LES is an unsteady solution technique that directly computes 
the large-scale turbulent structures in the flow. Because the large structures carry the vast majority of the 
turbulent kinetic energy, directly computing them should improve predictions. The effect of structures that 
are too small be be resolved by the computation is modeled or neglected. While LES has the potential 
for improved accuracy, up to this point RANS methods perform as well or better at a much lower cost for 
compressible jets. 6 A recent survey comparing several different LES methodologies for the prediction of jet 
flows was done by Bodony and Lele. 7 In general the methods do a poor job of predicting the end of the 
potential core, centerline velocity profiles and centerline turbulence intensities. 

The key to an accurate LES is the proper resolution of the large-scale turbulent structures which dominate 
the flow. The resolution of these structures is determined by a combination of the numerical scheme and 
computational grid used. Improved resolution can be achieved by adding additional grid points or increasing 
the resolution of the scheme. To quantify the ability of numerical schemes to resolve and convect wavelike 
structures many researchers have used a Fourier analysis based technique. 8 This technique has been a 
valuable tool in the development and characterization of numerical schemes. Recent work has attempted to 
extend this work to jet LES. 9 


NASA/TM— 201 1-216833 


2 



High resolution schemes are defined here as having the ability to resolve structures containing large wave 
numbers with fewer grid points than the typically used second- and third-order schemes. High resolution 
schemes may, but do not necessarily, have a high order of accuracy. Commonly used high resolution schemes 
are compact schemes 10 and Dispersion Relation Preserving (DRP) schemes. 11, 12 These schemes are more 
efficient for turbulent simulations; they can achieve a given level resolution with less computer resources than 
low-order schemes. A drawback of high resolution schemes is the increased complexity required in the flow 
solver. Large difference stencil sizes make it difficult to implement the numerical scheme near computational 
boundaries and more information must be passed between grid interfaces. This makes using high resolution 
schemes difficult for realistic problems with complex shapes. 

This paper describes a new LES code and its application to jet flows. The code incorporates high 
resolution numerical schemes for efficient simulations of turbulent flows. Details on the numerical schemes, 
boundary conditions and parallel processing options are provided. Several features that facilitate the solution 
of flows around complex geometries, generalized curvilinear coordinates, internal grid boundaries, hole cutting 
and grid blocking, are included and their implementation with respect to the spatial differencing schemes 
is discussed. The code is applied to a Mach 0.9 turbulent compressible jet, a standard test case, and its 
solutions are compared against experimental data. Mean velocities and turbulence intensities are compared 
to the experiment with good results. 


II. Code Description 

The code used in this study, WRLES (Wave Resolving Large-Eddy Simulation), is a special purpose large- 
eddy simulation code that uses high-resolution temporal and spatial discretization schemes to accurately 
simulate the convection of turbulent structures. The code solves the compressible Favre- filtered Navier- 
Stokes equations. The code is written entirely in Fortran 90 and utilizes both Message Passing Interface 
(MPI) libraries 13 and OpenMP compiler directives. 14 It was developed at the NASA Glenn Research Center 
for the study of turbulent jets, but can be applied to other flows. 


A. Governing Equations 

The basis of large-eddy simulation is the separation of the large and small scale turbulent fluctuations. The 
large-scale turbulence is computed directly using the Navier-Stokes equations. The small scale turbulence is 
modeled. Since the large-scales carry the vast majority of the turbulent kinetic energy the technique offers 
promise for accurate turbulent simulations. Since the small scale turbulence serves to dissipate the large 
scales and is isotropic, a simple model can be constructed to impose the effect of the small scales on the rest 
of the flow. 

A filtering process is applied to the continuity, momentum, and energy equations. The resulting equa- 
tions are comprised of resolved and unresolved terms. The resolved terms in the filtered equations directly 
correspond, in form, to the unfiltered equations. The resulting LES expressions for conservation of mass, 
momentum, and energy are: 
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The overbar, ( ) represents a spatially filtered quantity and the tilde, ( ) represents a Favre-averaged 
quantity; /= (pf)/p 

The unclosed terms from the LES momentum and energy equations are the sub-grid scale stress, r/? s , and 
the sub-grid scale heat flux, q S9S . In a traditional LES approach a sub-grid model would be used to compute 
r-- s and q S9S . The Smagorinsky sub-grid model is included in the code. 15 However, many practitioners 
forego sub-grid modeling, relying on the dissipation implicit in the numerical scheme to dissipate the energy 
at the small-scales. This approach is sometimes called Implicit LES (ILES). 
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B. Numerical Method 


The Favre filtered Navier- Stokes equations are discretized and solved using the explicit numerical methods 
described below. 


1. Time Discretization 

The code uses a family of explicit Runge-Kutta time stepping schemes written in a general M - stage 2 -N 
storage formulation. 16 

dfm == Qtmdfm— 1 T AtD (/ m _ i) (4a) 

fm = fm—1 + Pmdfm (4b) 

for m = 1 ... M, and where the initial value at the start of the iteration is /o = f n and the final value is 
/m = / n+1 . The coefficient aq is typically set to zero for the algorithm to be self-starting. 

The operator D is the spatial finite difference operator. One- stage is required for each order of accuracy 
desired. Additional stages can be used to increase accuracy or alternatively, to reduce dispersion error. The 
temporal discretization can be varied by simply changing the number of stages and the coefficients, a m and 
Pm- Table 1 lists the schemes that are currently implemented in the code. 


Scheme 
Euler forward 
Second order 
Williamson 16 
Gottlieb & Chu, TVD 17 
Carpenter & Kennedy 18 
Carpenter & Kennedy 18 
Stanescu & Habashi 19 
Stanescu & Habashi 19 


Number of Stages 
1 
2 
3 

3 

4 

5 

5 

6 


Order of Accuracy 
1 
2 
3 
3 

3 

4 
4 
4 


Table 1. Low-dispersion Runge-Kutta schemes implemented in the code 


2. Spatial Discretization 


Central differencing is used for the spatial discretization because of its non- dissipative properties. This helps 
ensure the accurate convection of turbulent structures. The central difference stencil can be written for an 
arbitrary stencil size 

dj_ 

dx 


— ^ ^ ^%[ a n(fi+n fi—r 


( 5 ) 


where the half width of the stencil, the number of points on either side of the central point, is N . The total 
number of points in the stencil is 27V + 1. The spatial discretization can be varied by simply changing the 
width of the stencil and the coefficients, a n . Standard stencils from 2nd- to 12th-order are included. In 
addition Dispersion Relation Preserving (DRP) stencils from Tam, 11 7-point, and Bogey & Bailly, 12 7-, 9- 
and 13-point, are included in the code. 

Standard central difference stencils are designed to minimize the truncation error for a given stencil 
size. In DRP schemes, the coefficients in the stencil are chosen to minimize the dispersion error of the 
scheme, rather than the truncation error. DRP schemes are designed to accurately capture and convect 
large structures, which are well resolved by the grid, and are generally a better choice for LES. 


3. Filtering 

The lack of dissipation in central differencing makes the schemes unstable. To ensure a stable solution 
without adversely affecting the resolving properties of the scheme, solution filtering is used. This is a 
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low-pass filter that leaves the low-wavenumber well-resolved structures untouched and removes the high- 
wavenumber unresolved structures that can cause instability. The filter must be properly matched to the 
differencing scheme, so that the filter removes only those waves that are not properly computed. 

N 

fi= fi+o- ^ b n fi +n (6) 

n=-N 

For the standard central difference schemes the filters of Kennedy and Carpenter 20 are implemented in the 
code. Bogey and Bailly developed filters to match their DRP stencils and are included. 


4- Properties of the Numerical Scheme 

Fourier analysis has been used to evaluate the ability of numerical schemes to resolve wave motion. 8 The 
technique analyzes the numerical scheme on the one-dimensional convection equation and can provide the 
scheme’s behavior in terms of dissipative and dispersive errors. The effect of the temporal discretization is 
not considered here. 

Fourier analysis shows that central difference schemes have no inherent dissipation, and thus are ideally 
suited for LES calculations. They do however produce dispersive errors. Figure 1(a) shows the phase error, 
written as the ratio of numerical phase speed to actual phase speed, c*/c, versus wavenumber per grid 
spacing, (kAx), for both standard schemes and the DRP schemes. If one selects a maximum acceptable 
phase error for a simulation then the maximum resolvable wave number per grid spacing, (^Ax) maa ,, is 
determined for each stencil. The maximum acceptable error used was 5 • 10 -4 . 

The damping function of the filter can also be expressed in terms of wave number per grid spacing. The 
damping functions, D(kAx), are shown in figure 1(b) for the standard filters of Kennedy and Carpenter 20 
and the DRP filters derived by Bogey and Bailly to match their difference stencils. The damping is zero 
until it reaches a “cutoff” wave number, where it rapidly increases, effectively removing all the structures of 
higher waves numbers. Since filters typically also utilize a coefficient, cr, to regulate the amount of damping, 
the total amount of dissipation is (jD{k Ax). The cutoff wave number, {nAx) cut is the wave number where 
the total amount of dissipation exceeds 5 • 10 -4 , the same error level used for the differencing schemes. A 
typical damping coefficient of 0.2 is used here. 




(a) Error in phase speed for the finite difference stencils 


(b) Damping functions of the explict filters 


Figure 1. Properties of the numerical scheme 


The maximum resolvable wave number for the 8th, 10th, and 12th order standard schemes and the 9-, 
11- and 13-point DRP schemes is plotted against the stencil size as solid lines in figure 2. The cutoff wave 
number of the corresponding filter is plotted as a dashed line. For a stable solution, the cutoff wave number 
of the filter should be less than the maximum resolvable wave number of the numerical scheme. This ensures 
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that all of the structures within the domain are properly resolved. For the standard schemes, the filter’s 
cutoff is too high, leaving improperly resolved structures in the flowfield. For these schemes, choosing the 
next lower order filter effectively shifts the filter curve to the right and should insure a stable solution. Figure 
2 clearly illustrates superior resolution offered by the DRP schemes. 



Figure 2. Cutoff wave numbers for schemes and filters 


C. Geometry Handling 

Many high resolution/high-order of accuracy CFD solvers are limited to canonical problems with very simple 
geometries. The WRLES code contains several features that allow the solution of flows about moderately 
complex geometries. 

1. Generalized Curvilinear Coordinates 

The governing equations are cast in generalized curvilinear coordinates allowing the computational grids to 
stretch and deform, fitting the geometric surface. The code has options to solve these equations in either the 
strong conservation form or the chain rule form. 21 The strong conservation form is used in most codes. 22 
However, Hixon et al. 23 recently showed that in practice, the chain rule form of the equations is more 
accurate than the strong conservation form when the metric terms are computed numerically. There are 
two practical considerations when choosing the form of the equations to solve; 1) the strong form is more 
computationally efficient, but 2) the chain rule form is more stable when using the centerline differencing 
condition for o-grid singularities. Grid metrics for the coordinate transformation are computed using the 
same differencing scheme that is selected for the flow solver. 

2. Internal Boundaries and Holes 

To allow greater flexibility in generating grids, boundary conditions can be specified on any portion of 
any grid surface in the domain. In addition, holes within the computational domain can be created. The 
combination of these features allows for complex geometric features to be embedded within a grid block. 

3. Near Boundary Stencils 

As you approach the boundaries, the full width of the difference stencil can not be maintained. To preserve the 
non- dissipative and non-dispersive characteristics of the scheme, the stencil width is reduced in the direction 
of the boundary and central differencing stencils are maintained. This approach reduces the resolution of the 
scheme near the boundaries and care must be taken to not adversely affect the solution. Wherever possible 
the boundaries should be placed as far from the region of interest as possible. For the jet problem, only the 
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very start of the mixing layer will be affected by the smaller stencils. Figure 3 illustrates the reduction of 
the stencil half widths, TV, near an internal boundary representing and embedded object. The boundary is 
represented by the red line and the hole points are shown by the x’s. 



(a) i-direction (b) j-direction 


Figure 3. Stencil half widths, N , near an embedded object 


4- Grid Blocking 

The code is also capable of computing grids made up of multiple grid blocks. In order to maintain consistent 
numerical resolution across the block boundaries, the blocks must be point to point matched and share 
several common planes. The number of common planes required is 27V, one less than the stencil width. This 
restriction ensures that every grid point in the domain is computed using a full central difference stencil. 
Block interfaces for stencil widths of 5 and 13 are currently included in the code. A schematic showing the 
overlap requirement and the resulting stencil half widths at each point for a 13-point stencil is given in figure 
4. 



Figure 4. Stencil half widths, N, at overlapping block boundary 


D. Boundary Conditions 

Boundary conditions are available to simulate a wide variety of flows. Because the width of the central 
difference stencil is reduced as the boundary is approached, reducing the resolution, most of the boundary 
conditions do not attempt to maintain high resolution. This is true for wall, inflow and outflow boundaries. 
High resolution is maintained for periodic, cylindrical centerline and block interface boundaries. The formu- 
lation of the boundary conditions is standard within the aerospace community. Table 2 lists the available 
boundary conditions. 
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Boundary Condition Condition(s) enforced 


subsonic inflow 
supersonic inflow 
characteristic boundary 
subsonic outflow 
supersonic outflow 
inviscid (slip) wall 
viscous (no-slip) wall 
singular axis 
periodic 
block interface 
centerline differencing 
exit zone 


Po, T 0 

P, pu, pv, pw , pe t 

Moo, a 

Poo 

none 

dui/dn = 0 

Ui = 0, dq/dn = 0 


Table 2. Standard boundary conditions 


An o-grid that forms a cylindrical coordinate system is the natural choice for a round jet problem. 
However this grid type creates a singularity on the centerline causing two problems; 1) a boundary condition 
must be applied in the center of the domain where the flow should be computed and 2) the small grid cells at 
the singularity create numerical instabilities. The centerline differencing condition removes both problems 
by removing the singularity, increasing the cell size, and creating a difference stencil across the centerline so 
that the equations can be solved there. The formulation is taken from Hixon 24 and the concept is illustrated 
in figure 5. The grid is constructed with a cylindrical void on the centerline to avoid a collapsed surface. 
The cell spacing is set such that the diameter of the void is equal to the adjacent radial grid spacing. Then 
a stencil is constructed in such a manner that the points on the opposite side of the singularity are used. 



Figure 5. Sample differencing stencil across the centerline 


The exit zone condition applies additional damping near an outflow boundary to remove outgoing waves 
that would otherwise reflect and contaminate the solution. The damping is introduced through grid stretching 
and additional low-order filtering whose strength is gradually increased as the boundary is approached. 
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E. Parallel Processing 

WRLES is written using hybrid parallel constructs to take advantage of multi-processor and multi-core Linux 
cluster systems. The hybrid parallel implementation uses both Message Passing Interface (MPI) libraries and 
OpenMP compiler directives. The grid can be divided into multiple blocks and computation of each block 
can be distributed to the nodes of the Linux cluster using the MPI libraries. OpenMP compiler directives 
parallelize the computation at the loop level across the multiple processors and multiple cores within each 
computer node. The requirement of grid overlap adds additional grid points to the computation when a 
domain is subdivided for parallelization. This limits the parallel scalability for large numbers of processors. 

III. Simulation of a Mach 0.9 Jet 

The WRLES code was used to obtain flowfields for a Mach 0.9 turbulent cold jet, a standard test case 
for jet LES. The configuration used for this study was the 2 inch diameter Acoustic Reference Nozzle (ARN) 
tested at the NASA Glenn Research Center by Bridges and Wernet. 25 The experimental data was obtained 
using Particle Image Velocimetry (PIV). The data was obtained on a streamwise plane through the jet 
centerline downstream of the nozzle exit from 1 < x/Dj > 25 and —1.4 < r/Dj > 1.4. Three components of 
velocity and turbulence intensities are used for comparison to the calculations. The jet Reynolds number, 
R e j = pjUjDj / fij, from the experiment was approximately 2 million. In order to reduce the range of 
turbulent scales in the simulation, the Reynolds number of the simulation was reduced to 50,000. 

A. Computational Approach 

1. Computational Grid 

A blocked structured grid was created using the Gridgen software package. 26 The computational grid (figure 
6) models the internal and external nozzle geometry from plenum to exit. The majority of grid points are 
placed in the area of interest, downstream of the nozzle exit in the plume region. The grid extends 40 jet 
diameters, Dj , downstream of the nozzle exit and 30 diameters radially from the centerline. The nozzle lip 
thickness was 0.015ZL, (0.03 inches) and 30 grid points were used to model it. For the viscous walls, the grid 
spacing at the wall was set to 3 • 10~ 4 Dj (6 • 10 -4 inches). 



Figure 6. Computational grid, xy-plane (every other grid line removed) 


The grid is constructed in 3 grid blocks; interior nozzle flowpath, exterior nozzle flowpath and nozzle 
plume. To ensure that each grid point is computed using the maximum stencil width possible the nozzle 
plume grid block is overlapped into the interior and exterior nozzle blocks. The grid points that fall in 
between the two upstream blocks are designated as hole points and not computed. Figure 7 shows a detail of 
this region. The overlapped region, hole points (in purple) and coincident points can be seen in the exploded 
view. 
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x/D. 


(a) 3-block grid around nozzle lip 


(b) exploded view, showing overlap and i-blanked region 


Figure 7. Modeling of nozzle lip 


Three different grid densities were examined (table 3). The initial coarse grid was constructed based on 
best practices for a Reynolds Averaged Navier-Stokes (RANS) simulation. For this grid, the smallest grid 
spacings were in the radial direction, due to the clustering near solid walls, and the largest were in the axial 
direction. The ability to resolve a three-dimensional turbulent eddy is limited by the largest grid spacing in 
each of the three directions. Therefore the first refined grid maintained the radial grid spacing and added 
points in the axial and azimuthal directions. The second refinement added points only in the axial direction. 


Grid 

Interior of Nozzle 

Exterior of Nozzle 

Nozzle Plume 

Total 

coarse 

45 x 55 x 102 

45 x 65 x 102 

196 x 148 x 102 

3,509,616 

medium 

45 x 55 x 132 

45 x 65 x 132 

294 x 148 x 132 

6,456,384 

fine 

45 x 55 x 132 

45 x 65 x 132 

435 x 148 x 132 

9,210,960 


Table 3. Dimensions of the computational grid 


2. Boundary Conditions 

Total pressure and total temperature, consistent with a Mach 0.9 cold jet, were specified at the inflow of the 
nozzle plenum. Freestream conditions were specified using a farfield characteristic boundary condition. The 
experiment had no freestream velocity, but a very low Mach number, 0.05, was specified at the computation’s 
farfield boundary to maintain stability. Static pressure was specified at the outflow and an exit zone was 
initiated at x/Dj > 30; the grid spacing is gradually increased and a sixth-order filter is applied to remove 
spurious reflections. The solution was computed at the jet centerline using the singularity differencing 
boundary condition. 

3. Solution Parameters 

Temporal discretization was done using the four-stage third-order low-dispersion Runge-Kutta scheme de- 
veloped by Carpenter and Kennedy. 18 Spatial discretization was done using the 13-point DRP scheme 
developed by Bogey and Bailly. 12 The filter corresponding to this scheme was used for damping with a 
coefficient of 0.5. 
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For these calculations the implicit large-eddy simulation (ILES) approach was used. In this approach, 
no model is used for the sub-grid scale turbulent stresses. The solution relies on the dissipation of the filter 
to represent the sub-grid scales. 

The solution does not utilize any artificial means to generate or initialize turbulence in the flowfield. 
The grid for the nozzle boundary layers is not resolved enough to sustain turbulent structures, resulting in 
laminar flow at the nozzle exit. At the nozzle exit, vortex shedding from the nozzle lip introduces large-scale 
unsteadiness into the jet mixing layer, which rapidly transitions into turbulence. 

The initial 3 block grid was subdivided into 16 blocks for parallel processing. The code was run on a 
linux cluster using one master node and 16 dual-processor dual-core worker nodes. The block interface data 
was shared between nodes via MPI libraries over an Infiniband network. The computational work on each 
node was parallelized across the four processor cores using OpenMP. 

4- Time Advancement 

The solution was run by specifying a CFL number. The time step at each grid point is computed using 
that CFL, and the minimum value for the entire domain is used to advance the solution everywhere. The 
maximum stable CFL for this calculation was 0.7, which corresponds to a physical time step of approximately 
2.45 • 10 -8 seconds. The solution was saved at a Nyquist rate of 2.5 • 10 -5 seconds, yielding a maximum 
frequency of 20kHz. A total of 2,048 solutions, representing 0.0512 seconds of physical time (7.7 flow through 
times) were saved for post processing. 

B. Results 

Instantaneous Mach number contours for the medium grid are shown in figure 8. Upstream of the nozzle 
exit the flow is steady. Downstream, the flow is dominated by turbulent mixing. Immediately downstream 
of the nozzle lip relatively large structures are seen; the structures are on the order of the the nozzle lip 
thickness. In general the size of the structures increases with distance from the nozzle exit and grow to be 
on the order of the nozzle exit diameter. The thin void along the centerline of the domain is an artifact of 
the centerline boundary condition. 



0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 


Figure 8. Instantaneous Mach contours, M 


1. Prediction of the Mean Flowfield 

The prediction of the mean velocity along the jet centerline for all three grids is compared to experiment in 
figure 9. Overall, agreement with experiment is very good for all three grid levels. The data shows that the 
predicted end of the potential core (the region of jet flow unaffected by viscous effects) shifts downstream 
with the initial grid refinement, more closely matching the data. Further refinement shifts the end of the 
potential core slightly back downstream. Downstream of the potential core, the decay rate of the centerline 
velocity increases with grid refinement. This is most likely due to the additional turbulent structures that 
are resolved by the finer grids. 
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Figure 9. Centerline velocity profile 


Radial profiles of mean velocity at four axial locations, x/Dj = 1, 4, 7 and 10, are shown in figure 10. 
Again, agreement with experiment is very good. The profiles indicate that the computations predict a thicker 
initial shear layer thickness than the experiment. Grid refinement appears to reduce the thickness, moving 
the predictions closer to the data. 


Coarse Grid 


Medium Grid 


Fine Grid 


□ Experiment 



Figure 10. Radial profiles of axial velocity 


To help quantify jet mixing a shear layer thickness parameter was defined. The inner edge of the shear 
layer was defined as the point where the velocity was 90 percent of the jet velocity and the outer edge was 
defined as the point where the velocity was 10 percent of the jet velocity. 


b = r 


— r 


U = 0.1Un 


U = 0.9Un 


( 7 ) 


The thickness, 6, based on the edge points was computed from the jet exit to the end of the potential core 
and is shown in figure 11. The experiment indicates a near linear growth rate from the jet exit. All three 
computations show that the shear layer grows rapidly for the first diameter and then relaxes to a growth 
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rate similar to the experiment. The initial thickness does not improve with grid resolution. This could 
be caused by; 1) the the lack of turbulence entering the shear layer from the nozzle boundary layers, 2) 
the fact that even the finest grid can not resolve the small scale turbulence here, or likely a combination 
of the two. Downstream the growth rates improve with grid resolution and the fine grid solution growth 
rate approaches the experimental data. For the coarse grid, the error in the initial shear layer thickness 
and the slow downstream growth rate compensate for each other. This is the reason the mean centerline 
velocity profile appears correct for this case. For the fine grid the initial shear layer is still too thick, but the 
downstream growth rate matches the experimental growth rate closely. The poor prediction of the initial 
thickness is the reason the fine grid centerline velocity decays too early. 



Figure 11. Shear layer thickness 


2. Prediction of Turbulent Statistics 

Contours of the three components of turbulence intensity are shown in figures 12 - 14. The three levels of 
grid resolution are compared to experimental data. For all components of intensity and all levels of grid 
resolution the computations show a high level of turbulence just downstream of the nozzle lip and persisting 
for 2 nozzle diameters. The experimentalists were unable to take accurate measurements in the first jet 
diameter downstream of the nozzle exit. Nevertheless based on the existing data one can infer that the high 
predicted levels of turbulence near the nozzle lip are nonphysical. Because the nozzle boundary layers are 
not properly resolved, large laminar-like vortices are shed. The situation is exacerbated by the lack of grid 
resolution in this region necessary to resolve small turbulent structures. These large vortices serve to disturb 
the flow and trigger the flow to transition to a realistic turbulent state. 

The contours show that the predictions of turbulence intensity are greatly improved with increasing grid 
resolution. For all three components of turbulence intensity the fine grid predictions beyond x/Dj = 5 are 
very good. 

The discrepancy near the nozzle lip can also be seen on line plots of turbulence intensity along the jet lip 
line, r/Dj = 0.5 (figure 15). A sharp peak of about twice the expected intensity in the downstream region is 
present near x/Dj = 0.5. For the coarse grid, levels of axial intensity are overpredicted and the radial and 
azimuthal values are underpredicted. Grid refinement reduces axial levels and increases radial and azimuthal 
levels improving agreement with experiment for all components. 

Axial, u'/Uj, and radial, v'/Uj , turbulence intensities on the jet centerline are shown in figure 16. On the 
baseline grid, the axial intensity levels are over predicted and the radial intensity levels are under predicted. 
Grid refinement improves both predictions, decreasing the axial intensity and increasing the radial intensity. 
Even though the flow at the centerline upstream of the end of the potential core is inviscid, the effect of 
the large structures at the nozzle lip are felt on the centerline as evidenced by a peak near x/Dj = 0. The 
disturbances from the jet lip are propagated to the centerline through pressure fluctuations. 


NASA/TM— 201 1-216833 


13 



(b) medium grid 



x/Dj 
(c) fine grid 


1.5 

1 

0.5 

0 



0 


5 


10 


15 


x/Dj 


(d) experiment 








0 

0.04 

0.08 

0.12 

0.16 

0.2 


Figure 12. Contours of axial turbulence intensity, u' /Uj 
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Figure 13. Contours of radial turbulence intensity, v'/Uj 
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Figure 14. Contours of azimuthal turbulence intensity, w' /Uj 
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Figure 15. Lip-line turbulence intensities 
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Figure 16. Centerline turbulence intensities 


Radial profiles of turbulence intensity at x/Dj = 1, 4, 7 and 10 are shown in figure 17. The trends seen 
in figures 12 - 16 are further quantified here: 1) the LES overpredicts the turbulence levels near the nozzle 
lip, 2) further downstream the turbulence levels are slightly underpredicted and 3) the predictions improve 
universally with grid refinement. 

The turbulence intensity data clearly show that grid refinement greatly improves the LES predictions. 
This fact was not evident in the mean velocity data, and the opposite conclusion could have been drawn. 
The effect of two competing sources of error, the poor resolution of the initial shear layer and the predicted 
growth rate downstream compensated for each other in the mean flow data. 

It is also important to note that the prediction of all components of turbulence intensity were improved by 
refining the axial and azimuthal grid spacing for the first refinement and only the axial spacing for the second 
refinement. In other words, improving the grid resolution in the coarsest direction has a significant effect. 
This illustrates the importance of resolving the turbulent structures in all three computational directions 
and suggests that isotropic grid cells are ideal for LES. This is a significant shift from the practices of RANS 
analyses where grid refinement is done to resolve mean flow gradients. 

3. Prediction of Turbulent Spectra 

Turbulent spectra were computed at three locations in the jet mixing layer. The points were located on the jet 
lip line, r/Dj = 0.5 at three axial locations, x/Dj = 4, 7, and 10. At each of these points, the three velocity 
components for each saved solution were extracted, yielding a time history of the velocity. The data was 
converted to instantaneous turbulent kinetic energy and was processed using a Fast Fourier Transform. This 
yielded a turbulent spectrum in the frequency domain. The flow was assumed to be statistically stationary 
and Taylor’s hypothesis 27 was used to convert the frequency domain to the spatial domain. 

Plots of the turbulent spectra for both grids are shown in figure 18. The majority of the turbulent energy 
is contained in the small wave number /large- scale structures. A dashed line with a —5/3 slope indicating 
the inertial sub-range is fitted to each spectra. The inertial sub-range is not well resolved (some may argue 
whether or not it is present in the simulation at all). Beyond this point the energy decays very rapidly, most 
likely the effect of the filter. The spectra show that the energy in the larger scales increases with downstream 
location. This is consistent with the fact that the turbulent structures grow in size as they move downstream. 

Grid refinement increases the energy in the small scales and shifts the start of the inertial sub-range, the 
—5/3 slope, to higher wave numbers. The apparent extent of the inertial sub-range increases slightly with 
grid refinement. It is interesting to note that the turbulence field is well predicted on the fine grid despite a 
lack of resolution of the turbulent inertial sub-range. 
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Figure 17. Radial profiles of turbulence intensity 
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(c) x/Dj = 10 


Figure 18. Turbulent kinetic energy spectra 
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IV. Summary and Conclusions 


A code for large-eddy simulation of turbulent jet flows was developed. The code employs low-dispersion 
Runge-Kutta time discretization and central differencing spatial differencing. It is coded in a general manner 
such that any number of temporal and spatial schemes that can be written in the given form can be used. 
Fourier analysis of the spatial discretization scheme was used to characterize the ability of the scheme to 
resolve waves on a computational grid and it was shown that the Dispersion Relation Preserving schemes 
are superior to standard central differencing schemes for the resolution of turbulent structures. Details on 
the implementation of the numerical scheme for generalized curvilinear coordinates, blocked structured grids 
and internal boundaries are provided. 

The code was applied to the flowfield of a Mach 0.9 cold jet. A four-stage third-order low-dispersion 
Runge-Kutta scheme was used for the temporal discretization and the 13-point Dispersion Relation Preserv- 
ing scheme of Bogey and Bailley 12 was used for the spatial discretization. A filter designed to match the 
spatial discretization scheme was used to maintain stability and to represent the sub-grid scale dissipation. 
Three different grids of 3.5, 6.5 and 9.2 million points were examined. The results of the computation were 
compared to the experimental data of Bridges and Wernet. 25 Excellent agreement was obtained for the mean 
velocities at all grid resolutions. Good agreement was achieved for the turbulence intensities for the baseline 
grid and the agreement improved significantly with grid refinement. The results represent an improvement 
over the previously cited RANS 2 and LES 7 results. The improvement with selected grid refinement suggests 
that grid cells that tend toward isotropic are ideal for predicting turbulent quantities. Results indicated that 
the prediction of the initial portion of the shear layer 0 < x/Dj < 2 is a limiting factor in the accuracy of 
the simulation and further research is needed to improve predictions in this area. 

Turbulent energy spectra were obtained at three points in the flowfield for both grids. The spectra 
appear reasonable, but lack a well defined inertial sub-range. A —5/3 slope turbulent decay can be seen for a 
small range of wave numbers. Despite a lack of resolution of the turbulent inertial sub-range, the turbulence 
intensities were well represented on the fine grid. 
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